TLDR: These updates resolve some issues with the quality control and vertical interpolation algorithms used for init_type="metgrid". There is not an obvious "best" method for quality control and processing of the origin model data.

One challenge arises from origin model points that are in close proximity, which can cause undesirable behavior in the final interpolation. WRF's real.exe tosses out the higher elevation point if two points from the origin model are within 500 Pa, which occurs often due to the surface and an output level having similar pressures. Here, we adopt a similar approach, however, since we are operating in a height-space instead of WRF's pressure-space, we assume a baroclinic atmosphere and roughly estimate the pressure delta. The threshold pressure delta for quality control is an input parameter (metgrid_proximity) defaulting to 500 Pa. Note that this ERF default value changed in this PR to match what is default for real.exe.

Most of the quality control difficulty arises due to potential differences between the origin model topography and that ingested by WPS's geogrid.exe. Even with the "high resolution" 30 second WPS dataset, there are occasional differences of hundreds of meters when processing GFS. Moving to SRTM or another high resolution DEM would likely result in more grid points noticeably deviating from the origin model's DEM. We want and need higher resolution topography, so what options do we have to handle this mismatch?

1. Use the profile as-is. A major problem here is that cases where the WPS DEM is significantly higher than the origin model DEM would result in unphysical, high near-surface velocities since they are sampled from a much higher altitude in the origin model. This should not be an issue in the initial conditions, but forcing the lateral boundaries with these high velocities near-surface could result in some unexpected behavior for otherwise innocuous situations. This is the default behavior of real.exe.
2. Use the profile as-is except we toss out any origin model data that is below the first 6 ERF vertical levels. This is the default behavior for ERF in this PR. The number of vertical levels considered here is an input parameter metgrid_force_sfc_k. The idea here is to emphasize the surface value and ensure it is used in the interpolation for at least the first 6 ERF vertical levels. This is what WRF's real.exe claims to be doing, but looking at output it seems that there is a bug and the profile is used as-is.
3. Set the lowest ERF vertical level using the "surface" values from the origin model. This may get around some of the problem described above, but now some vertical gradients may become unrealistic and problematic. This behavior can be enabled by setting metgrid_retain_sfc = true.
4. Manipulate the profile in other more invasive ways. Squish, stretch, crop, etc, with the intention of capturing the spirit of the near-surface region while not messing with the conditions aloft. This would be an artistic set of choices and would have its own challenges.

Above, the terrain height from the origin model (GFS) is nearly 1 km lower in elevation than the terrain height used in WPS. This mismatch results in higher velocities and lower vapor mixing ratio in the near-surface relative to the origin model.

Above, the terrain height from the origin model (GFS) is several hundred meters above the terrain height used in WPS. This mismatch results in a zigzag behavior near-surface in the vapor mixing ratio.

Our next challenge relates to origin model data that is below the origin model terrain, which is common since many NWP data sources are provided on pressure levels. This data is typically filled to provide some agreement with the above-surface profile, however, the fill approach is not clear and may vary by origin model. For GFS, it appears that the below-surface velocities are extrapolated from the first two levels above-surface; the potential temperature is constant and equal to that at the first level above-surface; and the vapor mixing ratio has a constant rate of change of 1.5 g/kg/km extending downwards from the first level above-surface. Do we trust this data? real.exe does. I'm skeptical. The origin model surface level often deviates from the trend of levels aloft, which does seem physical. Levels below the surface tend to agree with levels aloft, meaning that when interpolating with levels above AND below the surface, the surface value can throw a stick in the spokes so-to-speak and result in awkward zigzags. As of this PR, the ERF default is to use both the surface value (metgrid_use_sfc = true) and use levels below the origin model surface (metgrid_use_below_sfc = true), which is the same default behavior as real.exe. In my opinion, metgrid_use_below_sfc = false is worth considering as a default.

Here, the vapor mixing ratio and potential temperature display a zigzag behavior that arises when both the surface and points below the surface are used (metgrid_use_sfc = true and metgrid_use_below_sfc = true) and the origin model terrain height is lower than the WPS terrain height.

Same as above, however, the ERF profile is constructed using the surface point (metgrid_use_sfc = true) but not the points below the surface (metgrid_use_below_sfc = false). Note the potential temperature profile near-surface is now neutral rather than being strongly stable.

Another parameter at hand is the order of the interpolation scheme; controlled by lagrange_order in the domains section of WRF's namelist.input, and the ERF input parameter metgrid_order. Everything shown above used 2nd order, which seems adequate when processing mesoscale NWP data, such as GFS or HRRR. There appears to be a bug in real.exe when using odd interpolation orders that produces some nonsense values near the surface. Near the data limits real.exe seems to be using out of bounds values. ERF gets around this by reducing to linear interpolation if the interpolation stencil reaches out of bounds. If necessary, this behavior could be improved in the future to gradually reduce the interpolation order when approaching data limits.

Third order interpolation with lagrange_order=3 in namelist.input and metgrid_order=3 in inputs. Note the anomalous results for the lowest WRF level.

Opinions, comments, and questions are all welcomed. Please share your insights and thoughts on how we can improve!